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1 Introduction 

The observation of Ultra High Energy Cosmic Rays (UHECR), begun in the 1960s, represents 
a unique window opened on the most energetic acceleration phenomena in the Universe. 
UHECR are observed at extremely high energies up to 3 -i- 5 x 10^*^ eV and the determination 
of their characteristics is of paramount importance in unveiling their possible astrophysical 
sources and/or acceleration processes. One of the key points of their study is related to the 
propagation of UHE particles in intergalactic space. The study presented here is mainly 
devoted to this analysis, outlining a novel computation scheme to treat the propagation of 
UHE particles. 

The propagation of UHECR from the source to the observer is mainly conditioned by the 
intervening astrophysical backgrounds, such as the Cosmic Microwave Background (CMB) 
and the Extragalactic Background Light (EBL). Experimental observations of UHECR should 
be always compared with theoretical expectations in order to firmly determine the nature of 
such fascinating particles and, maybe, their sources. 

Several propagation dependent features in the spectrum can be directly connected with 
the chemical composition of UHECR and/or to the distribution of their sources [1-3]. Among 
such features, particularly important is the Greisen-Zatsepin-Kuzmin (GZK) suppression of 
the spectrum [1], an abrupt suppression of the observed proton flux due to the interaction 
of UHE protons with the CMB radiation field. The GZK suppression, as follows from the 
original papers [1], refers to protons and is a consequence of the photo-pion production 
process {p + 'Jcmb — > p + vr) suffered by these particles interacting with the CMB radiation 
field. The energy position of the GZK cut-off, as well as the flux behavior in its proximity, 
can be predicted theoretically with extreme accuracy [4]. 

In the case of UHE nuclei the expected flux also shows a suppression at the highest 
energies due to the photo-disintegration process on the CMB and EBL fields, with the pro- 
duction of secondary (lighter) nuclei and nucleons: A + jcMB,EBL — > {A — nN) + nN, where 
A is the atomic mass number of the nucleus. The energy position of the suppression in the 
spectrum depends on the nuclear species, mainly on its atomic mass number A, and on the 
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details of the astrophysical backgrounds [5]. Particularly relevant is the EBL field, which 
fixes the energy of the onset of the flux suppression [5] . 

Another important quantity that, in principle, could affect the fiux behavior at the 
highest energies is the maximum energy provided by the sources Emax- In a typical scenario 
of rigidity dependent acceleration, the maximum acceleration energy of nuclei is proportional 
to the same quantity for protons through the nucleus atomic number (charge) Z, being 
^max = ZEmax- Therefore, for sufficiently low Emax the UHECR fiux steepening at the 
highest energies could be directly linked with the nucleus charge following, in this case, a 
picture analogous to the "knee" behavior observed in the case of galactic CR [6]. 

As a general remark it should be stressed that UHE protons propagation is affected 
only by the CMB field, since its density is almost three orders of magnitude larger than that 
of EBL [7, 8] . Therefore at all energies the proton energy losses on CMB largely dominate 
over those on EBL [5]. In the case of UHE nuclei the photo-disintegration process on EBL is 
relevant because, in the Lorentz factor range T < 2 x 10^ [5], it has no CMB counter-part and 
it changes substantially the expected fiuxes. Finally, the pair-production process of nuclei 
over the EBL field, as in the case of protons, is negligible because always dominated by the 
CMB radiation field. 

The EBL radiation field suffers of several uncertainties mainly connected with its cos- 
mological evolution while the CMB field is analytically known at any red-shift. This fact 
explains why the propagation features in the spectrum of protons, such as the GZK-cut off, 
are less affected by uncertainties with respect to the ones relative to nuclei. 

From 1960s a flattening has been observed in the UHECR spectrum at an energy around 
3 ^ 6 X 10^8 eV, which was called "the ankle". This feature may be explained in terms of 
the pair-production dip [2], that, like the GZK steepening, can be directly linked to the 
interaction of protons with the CMB radiation. The dip arises due to the process of pair 
production suffered by protons interacting with the CMB field p + ^CM _b — >■ p + e+ -|- e~ [2] . 
It is present only in the spectrum of protons at energies in the range 2 x 10^*^ 10^^ eV. 
The pair production process arises also in the propagation of nuclei, although it doesn't leave 
any feature in the expected spectra, through the reaction A + ^cmb A + e'^ + e~ and it 
involves only the CMB field, being the only EBL radiation relevant at the highest energies 
where the photo-disintegration process kicks in [5]. 

If nuclei dominate the UHECR spectrum the behavior of the observed fiux, the ankle, 
could have a different explanation that has been proposed by Hill and Schramm [9]. They 
interpreted the observed ankle in terms of a two-component model; the low energy component 
being either galactic or produced by the Local Supercluster. A similar model was later 
considered also in [10]. In this case the ankle energy region corresponds to the transition 
between two different components. 

From the experimental point of view the situation is still unclear. The HiRes experiment 
shows spectral features consistent with the proton GZK suppression and the pair-production 
dip [11]. Coherently with this picture, the chemical composition observed by HiRes is proton 
dominated at all energies E > 10^* eV [12]. Recent data from Telescope Array [13] appear to 
confirm this framework. The situation changes if the Auger results are taken into account. 
The Auger energy spectrum [14] shows with high statistical accuracy the two main spectral 
features: ankle and high energy suppression, but the corresponding energies are shifted with 
respect to the HiRes energies by about 25%. However this shift could originate from the 
different energy scales of the experiments whose systematic uncertainty is of the same size. 
The most discrepant outcome is in the possible interpretation of the mass composition from 
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the elongation rate data [15] which may be interpreted as a transition from hght to a heavier 
composition at energies E > A x 10^^ eV. This puzzling situation, with different experi- 
ments favoring different scenarios, shows the importance of a systematic study of UHECR 
propagation in astrophysical backgrounds. 

This paper describes a new Monte Carlo (MC) simulation code, SimProp^ , developed for 
the propagation of UHE particles (protons and nuclei) through astrophysical backgrounds. 
In designing such new scheme, which is not the first in this field of research [16-27] , we have 
focused on a tool which can provide a fast and reliable analysis of the predictions on the 
spectrum and chemical composition changing the background characteristics. 

In its current implementation SimProp uses a simplified nuclear model and a mono- 
dimensional treatment of the propagation, i.e. particles are propagated only in red-shift from 
the source to detection. More complete nuclear models and three-dimensional effects caused 
by the actual source distribution and the interaction of UHE particles with intergalactic 
and/or galactic magnetic fields will be included in further developments of the code. 

In the present paper we will present a systematic comparison of the SimProp results 
with the results of other MC schemes [16, 28] and with the analytical solution of the UHECR 
transport equations. The paper is organized as follows: in section 2 we discuss our theoretical 
treatment of the propagation of UHE particles in astrophysical backgrounds, in section 3 we 
introduce the layout of our MC code and its input-output, in sections 4 and 5 we compare 
the results of SimProp with other computation schemes and with the Auger observations on 
the spectrum respectively. Finally, conclusions and a discussion on future developments of 
the code take place in section 6. 

2 UHE Cosmic Ray Propagation 

The propagation of charged particles (protons or nuclei) with energies above 10^'^ eV through 
astrophysical backgrounds can be suitably studied taking into account the main channels of 
interaction that, as already anticipated in the introduction, are: 

• protons - UHE protons interact only with the CMB radiation field giving rise to the two 
processes of pair production and photo-pion production. We neglect their interaction 
on EBL as discussed in the introduction. 

• nuclei - UHE nuclei interact with the CMB and EBL radiation fields, suffering the 
process of pair production, in which only CMB is relevant, and photo-disintegration, 
that involves both backgrounds. While the first process conserves the nuclear species, 
the second produces a change in the nuclear species, extracting nucleons from the 
nucleus [5, 29]. 

In the energy range E ~ 10^^ -r- 10^^ eV the propagation of UHE particles is extended 
over cosmological distances with a typical path length of the order of Gpc. Therefore we 
should also take into account the adiabatic energy losses suffered by particles because of the 
cosmological expansion of the Universe. 

The computational scheme used to handle the propagation of charged particles in Sim- 
Prop is based on the kinetic approach proposed in [5]. The main ingredients of this method 

^The SimProp code here presented is available for the community upon request to: SimProp- 
dev@aquila.infn.it 
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are the continuous energy loss (CEL) approximation and the assumption of an exact conser- 
vation of the particle's Lorentz factor in the photo-disintegration process. Under the second 
hypothesis, namely neglecting the nucleus recoil in the interaction, we can easily separate the 
processes that change the Lorentz factor of the particle, leaving unchanged the particle type 
(pair and photo-pion production), from the processes that conserve it, changing the particle 
type (photo-disintegration) . 

The CEL approximation consists in assuming that particles lose energy (i.e. change 
their Lorentz factor) continuously. In the propagation through astrophysical backgrounds 
the interactions of UHE particles are naturally affected by fluctuations, with a non-zero 
probability for a particle to travel without losing energy. In the CEL approximation such 
fluctuations are neglected. 

In the case of proton propagation the CEL approximation has a negligible effect on the 
pair-production process, while in the case of photo-pion production it gives a deviation only 
at the highest energies {E > 10^^ eV) of the order of 10% with respect to the flux computed 
taking into account the intrinsic stochasticity of the process [2, 30]. Having this in mind, 
we have chosen in SimProp to handle nucleon propagation always under the CEL hypothesis 
using the analytic computation scheme presented in [2]. Moreover, in the case of nucleons we 
will not distinguish between protons and neutrons because [31]: (i) the photo-pion production 
process is an hadronic process that is essentially the same for protons and neutrons, (ii) the 
loss length of protons is always larger than the neutron decay length, apart from the extreme 
energies (fewxlO^*^ eV) where they become comparable [31]. Therefore in the following we 
will always refer only to protons. 

In the case of propagation of nuclei the energy losses due to the process of pair production 
can be simply related to the corresponding quantity for protons [5] 



with A being the atomic mass number of the nucleus and Z its atomic mass. Thus, using the 
results of [2, 30], we can use the CEL approximation also for the process of pair-production 
involving nuclei. 

The change in the Lorentz factor of the propagating particles is also linked to the 
cosmological evolution of the Universe. The expansion of the Universe causes an adiabatic 
energy loss to the propagating particles, that is (by definition) a continuous process common 
to protons and nuclei given by 



where H{z) = Hq^J (1 -|- z)^Om + is the Hubble parameter at redshift z in a standard 
cosmology with: Hq = 71 km/s/Mpc, = 0.24 and r^A = 0.72 according to WMAP data 



Let us now discuss the process of photo-disintegration of nuclei: this interaction changes 
the nucleus kind leaving its Lorentz factor unchanged. In the kinetic approach of [5] the 
process of photo-disintegration is treated as a decay process that simply depletes the flux of 
the nucleus A. Unlike the processes discussed so far, that are scarcely affected by fluctuations, 
the process of photo-disintegration could be much more affected by the stochasticity of the 
interaction. Therefore we have implemented our MC scheme only on this interaction process. 




,+ 



(2.1) 




(2.2) 



[32]. 
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which is simulated by computing the interaction time averaged over the density of the ambient 
photons: 



energy of the background photon in the rest frame of the particle, eo(^) the threshold of the 
considered reaction in the rest frame of the nucleus A, a the relative cross section, e the energy 
of the photon in the laboratory system and n-y(e) the density of the background photons per 
unit energy. Equation (2.3) is written using the Blumenthal approach [5, 33] and it refers to 
the specific reaction channel i, each characterized by a branching ratio, as reported in table 
1 and table 2 of [34]. The total inverse interaction time ta{T) can be obtained summing 
over the all possible photo-disintegration channels i. The photo-disintegration cross section 
as well as the relative branching ratios used in this work are taken from [34]. 

The dominant channels of photo-disintegration are single and double nucleon emission 
associated to the Giant Dipole Resonance (GDR) [34]. These processes are favored if the 
energy of the background photon in the rest frame of the nucleus is e < 30 MeV. At higher 
energies in the range 30 < e < 150 MeV a multi-nucleon emission regime takes over, while at 
energies e > 150 MeV the photo-disintegration cross section rapidly goes to zero [34]. 

Given the approximations described above, the SimProp computation scheme is a one 
dimensional algorithm in which only the red-shift z follows the " history" of the propagating 
particle. This approximation together with the Lorentz factor conservation in the photo- 
disintegration process justifies integrating over the photon density, as in equation (2.3), in- 
stead of generating the background photon parameters from their distribution. 

In our computation scheme the atomic mass number A uniquely tags the nucleus species. 
Following [34] we have chosen a list of nuclei from deuterium {A = 2) up to iron (A = 56) 
with one stable isotope for each atomic mass number A. This assumption is reasonable 
because, as discussed in [5, 17], the radioactive decay time to the line of stability is less 
than the one- nucleon emission photo-disintegration loss time^. In the case of mass values 
A = 54, 50, 48, 46, 40 and 36 there is more than one stable isotope but the laking of cross- 
section data for all stable isotopes makes impossible the computation, in these cases we 
have chosen the nucleus specie with a reasonable determination of the photo-disintegration 
cross section. A posteriori (see section 4), the agreement of our results with more refined 
computations schemes that take into account the effect of radio-active decays, such as the 
computations of [16], gives a solid justification to our approach. 

Following [5] we are not including the photo-pion production process for nuclei. This 
choice is motivated by the fact that nuclei photo-pion production is naturally suppressed 
because the energy of the photon in the nucleus rest system is A times lower than for a 
proton of the same energy. Differences with other simulation studies that include the photo- 
pion production for nuclei are observed only at Lorentz factors F > 10^^, as we will discuss 
in the next session, producing some differences over the corresponding energy in the spectra. 
However we note that at these high energies there is not enough statistics with current 
experimental data to compare different theoretical models and approximations. In future 
developments of the code we will come back to this approximation including also the tiny 
effects due to the photo-pion production process for nuclei. 

■^This statement is strictly correct for all nuclei species but three unstable nuclei: ^^Mn, ^^Al and ^^Be. 
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with A the atomic mass number and F the Lorentz factor of the interacting particle, e the 
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The average interaction time in Eq. (2.3) refers to the present epoch, for red-shift z = 0. 
Due to the expansion of the universe, both the background photon density and energy will 
evolve with red-shift. In the case of CMB this evolution is known analytically: the density 
changes as ucmb ^ (1 + z)^i^cmb and the energy as e — > (1 -|- z)e. The case of EBL is less 
clear: the EBL radiation is emitted by astrophysical objects at present and past cosmological 
epochs and subsequently is modified by red-shift and dilution due to the expansion of the 
Universe. The EBL energy spectrum is dominated by two peaks: one at the optical and the 
other at the infra-red energies, produced respectively by direct emission from stars and by 
thermal radiation from dust. At present there are only a few calculations of the EBL which 
include cosmological evolution, most notably [7] and [8]. In the present paper we mainly use 
the EBL as presented in [7] , which is a refinement of previous calculations [35] , based on the 
data from the Spitzer infrared observatory and the Hubble Space Telescope deep survey. In 
[7] the EBL photon density is found from 0.03 eV up to the Lyman limit 13.6 eV for different 
values of the red-shift up to z = 6, after which the EBL is supposed to be zero. 

3 Monte Carlo Layout 

The main ingredients to initiate the simulation process are: the red-shift of the source, the 
primary nucleus species and its injection energy at the source. The simulation code propa- 
gates particles in one dimension with only red-shift determining the particle evolution, since 
for a given cosmology there is a one-to-one correspondence between z and position of the 
particle. Consequently the particle is propagated in steps of red-shift. The MC follows the 
initial nucleus, secondary nuclei and protons produced at each photo-disintegration inter- 
action calculating their losses up to the observer, placed at red shift zero. As discussed in 
section 2, the propagation of UHE particles is based on the analytical scheme described in 
[2] for protons and in [5] for nuclei. The latter is modified for the use in a MC approach as 
described below. 

In this initial implementation the nuclear model adopted in SimProp is quite simple: 
following [34] we fix a list of nuclei that can be propagated, whose photo-disintegration cross- 
section is given in the same paper. Each nuclear species in the list is univocally identified by 
its atomic mass number A with steps of A A =1 starting from iron ^ = 56 down to beryllium 
A = 9. The unstable nuclei with 5 < A < 8 are excluded from the list and, for masses lower 
than A = 9, only helium A = 4, tritium ^4 = 3 and deuterium A = 2 are included. 

The values of the energy threshold for single or double nucleon emission are taken from 
[36]. In the case of isotopes with the same mass we choose the nucleus with the minimum 
value of the energy threshold for the emission of one nucleon (neutron or proton) . The energy 
threshold for the emission of two nucleons is chosen again as the minimum among the three 
different values for the emission of any pair of nucleons. 

Stochastic interactions must necessarily link a parent nucleus with one, or more, nuclei 
belonging to the list. The extracted nucleons (AA) are all treated as protons, as already 
discussed in the previous section. The choice of a nuclear model based on a limited list of 
nuclei relies on the hypothesis that a few representative processes can mimic a more complex 
description saving the computation time. The validity of this approximation, already used 
by different authors [16-18], will be verified a posteriori. 

The assumptions described above have an immediate consequence in the code layout, 
which is schematically sketched in figure 1: 
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Figure 1 . Flow chart of the simulation code SimProp. 



• nuclei follow a branch of the code where both continuous and stochastic processes 
occur. This is done in the method called PropagateNucleus through several steps each 
one determined by the actual occurrence of a stochastic process. Then a change in the 
nuclear species and the emission of protons occur. The steps are iterated up to the 
observer at z = 0. 

• Protons, which (within the CEL approximation) do not suffer stochastic interactions, 
are treated in the kinetic approach with a single step from their origin up to zero 
red-shift. This is performed in the method called PropagateProton. 

Let us now discuss the implementation of the stochastic treatment of the nuclei propa- 
gation. As described above, this is done in the method PropagateNucleus. The calculations 
of the energy evolution and of the photo-disintegration life-time (Eq. (2.3)) are performed 
step by step in red-shift. The survival probability as a function of red-shift and Lorentz 
factor of the nucleus A is: 

p(r,z) = exp(- r J; 

y ta{T,z) 

where z and z* are the values of the redshift of the current step (from z* to z). In the 



dt 
dz 



dz 



(3.1) 
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standard cosmology the term \dt/dz\ is given by 




Figure 2. Total photo-disintegration mean free path as a frmction of the Lorentz factor for iron 
nuclei with CMB and EBL backgrounds at red-shift z = 0. The black dashed line corresponds to 
the path length in the case of the EBL evolution of [7] (EBL 1), the red full line corresponds to the 
broken power law approximation of [7] (EBL 2), dots refer to the path length computed in [16]. 

In SimProp the intervals in red-shift which label the position of the particles have an 
exponentially decreasing size towards the generation point of the nucleus. This choice assures 
a higher accuracy in the evolution reconstruction near the production point. Given the 
probability (3.1) and the energy of the nucleus at a certain step in red-shift, the MC method 
is applied to select whether the nucleus interacts and to calculate the actual interaction 
point within the step. The MC method is also applied to decide the multiplicity of the 
photo-disintegration. To this end we have used the (average) cross sections with the different 
branching ratios as reported in tables 1 and 2 of [34]. 

The SimProp program is developed in C-|— 1-. The inputs needed by the code are: (i) 
the initial random seed; (ii) the number of events; (iii) the type of astrophysical background; 
(iv) the nucleus mass; (v) the minimum and maximum generation energy of the nucleus; (vi) 
the minimum and maximum generation red-shift of the nucleus. The simulation code can 
be run injecting at the sources either a fixed primary nucleus species or any distribution of 
nuclear masses. Details about the execution performances of SimProp are given in appendix 
A. 

The output of the simulation is stored in a ROOT [37] file recording the particles at 
each step of their propagation. The output is organized in branches containing the following 
information: (i) the branch of the propagation; (ii) the mass and the charge of the nucleus; 
(iii) the initial and final energy; (iv) the initial and final redshift; (v) the multiplicity of the 
interaction suffered by the current nucleus; (vi) the distance covered in the current step. 
The branch number zero refers to the primary nucleus. Nuclei and protons produced by 
photo-disintegration are traced branch by branch till they reach the observer at z = 0. 

The code is designed in such a way that any red-shift distribution of sources and any 
injection spectrum can be simulated. This is achieved drawing events from a flat distribution 
in the red-shift of the sources and of the logarithm of the injection energy. Once the event is 
recorded at z = the actual source/energy distribution is recovered through a proper weight 
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attributed to the event. As an example, let us discuss the case of uniformly distributed 
sources in co-moving coordinates with a power law injection spectrum. In this case events 
should be weighted with a factor 

Wz oc I , (3.3) 

with z the source red-shift. In the same way to generate a power law injection spectrum, 
with spectral index 7 and generation energy Eg, a weight 



WE cx El-"/ (3.4) 



has to be assigned to each event at z = 0. 
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Figure 3. Flux of Fe (lower curves) and all particle (upper curves) in the case of pure Fe injection 
with a power law index 7 = 2.2 and no energy cutoff at the source. Black: calculated with the EBL 
given in [7], EBL 1. Red: calculated with the analytical parametrization of the EBL background, 
EBL 2 (see text). 

Some of these inputs require a choice among options, depending on the specific needs 
of the user. As it is easy to understand, depending on the different options the performances 
of the simulation code could change. 

One of the input required by SimProp is connected with the EBL evolution model which 
is not analytically known, as discussed in section 2. The EBL evolution assumed in SimProp 
is the one given in [7] (EBL 1) or an analytical approximation of it based on a broken power 
law behavior (EBL 2); the latter choice assures a faster computation time. Other possible 
assumptions on the EBL evolution can be found in [8]. These different choices for the EBL 
evolution were tested with essentially the same final results: the only difference connected 
with this choice regards the execution time of the simulation. 

The total photo-disintegration path length for iron as function of the Lorentz factor is 
shown in figure 2. The black dashed line is calculated assuming the EBL evolution reported 
in [7] (EBL 1), while the red full line is obtained using its broken power law approximation 
(EBL 2) . The black dots represent the path length as computed in [16]. The differences 
between (EBL 1) and the simple analytical approximation (EBL 2) are limited to the low 
energy region, where the effect of photo-disintegration is not relevant being the corresponding 
path length of the same order of the Universe size (at Gpc scale). The differences with respect 
to the results reported in [16] are sizeable only at the highest energies where it is relevant 
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the effect of the photo-pion production and the weight of those nuclei that we have neglected 
in our simulation (see the discussion of section 2) . 

In any case what is important to discuss here is the effect on the observed flux on Earth 
due to the different possible assumptions on the EBL evolution. In figure 3 we plot the 
flux of primary iron expected on Earth in the case of an injection spectrum cx E~'^''^ with 
uniformly distributed sources in co- moving coordinates. Here and in the following all fluxes 
are normalized to the Auger spectrum [38] above 10^^'^ eV. The red triangles refer to the 
(EBL 1) while the black triangles to the simple broken power law behavior (EBL 2). Sizable 
differences in the iron spectra are recognizable only at energies below 10^9 eV. In the same 
flgure the all particle spectrum is shown, i.e. the sum of primary iron and all secondaries 
produced by photo-disintegration along the propagation path: in this case the effect of the 
different choices for the EBL is negligible. 




20 20.5 

logj^(E/eV) 



Figure 4. Flux of iron and secondary nuclei (A=50, 40, 30, 20, 10) at z = in the case of pure 
iron injection at the source with a power law injection index 7 = 2.2. Full squares correspond to the 
SimProp result while continuous lines correspond to the solution of the nuclei kinetic equation of [5]. 



4 Comparison with other propagation schemes 

In this section we discuss the comparison between the results of SimProp and other com- 
putations schemes based: (i) on a pure kinetic approach and (ii) on different MC schemes. 
In particular, since SimProp is based on the kinetic approach of [5], a comparison with the 
results obtained in such a scheme is of particular importance in order to assess the internal 
consistency of our MC code. To compare the SimProp results with other MC computa- 
tion schemes we have chosen the simulation by Allard et al. in [16] and the CRPropa [28] 
simulation code. Let us discuss separately the two cases. 

4.1 Kinetic Approach 

In this sub-section, the spectra obtained using SimProp have been compared with those 
calculated solving the kinetic equation associated to the propagation of nuclei [5] . To pursue 
such comparison, a pure iron injection with a power law injection of the type cx E^"' with 
7 = 2.2 have been assumed. The sources have been assumed to be homogeneously distributed 
in the red-shift range < z < 3. In flgure 4 the fluxes expected at z = are shown for 
iron and secondary nuclei produced in the photo-disintegration chain suffered by primary 
injected irons. The points refer to the SimProp results while the continuous lines to the 
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fluxes computed in the kinetic approach [5]. A good agreement between the two schemes 
is clearly visible in figure 4. At the highest energies the path- length of iron nuclei is very 
short (lower than few Mpc, see figure 2). Therefore, to achieve a good sampling in the 
MC simulation, higher statistics is needed; this is the reason for larger errors bars in the 
SimProp results at the highest energies and for their less good agreement with the solution 
of the kinetic equation. Notice also that the simulation used for this comparison has reduced 
statistics respect to the other figures and, more importantly, here secondaries are not grouped 
together. 




logj^(E/eV) 

Figure 5. Flux of secondary protons at z = in the case of pure iron injection at the source with 
a power law injection index 7 = 2.2. Full squares correspond to the SimProp result (instanteaneous 
photo-disintegration approach and standard SimProp approach) while the continuous line corresponds 
to the instanteaneous photo-disintegration applied to the kinetic approach of [5]. 

In the photo-disintegration chain of iron, among secondary particles, protons are also 
produced. As discussed in [5], the flux of secondary protons can be easily computed assuming 
an instantaneous photo-disintegration of the primary injected nucleus. In this case an iron 
nucleus once injected at the source with energy Eg is immediately destroyed into A = 56 
nucleons each of energy Eg/A. At large Lorentz factors this assumption is well justifled 
because the nucleus lifetime (2.3) is much shorter than all other relevant time scales of the 
problem (see also figure 2). 

In figure 5 we show the fiux of secondary protons expected at z = computed in a 
full SimProp simulation and assuming an instantaneous photo-disintegration of primaries. In 
figure 5 we have also computed the fiux of protons obtained in SimProp forcing to zero the 
nuclei path-length, to mimic the physics of the instantaneous photo-disintegration. Figure 5 
shows the expected behavior: the flux obtained by a full SimProp simulation is bounded from 
above by the flux obtained assuming an instantaneous photo-disintegration that coincides 
with the SimProp flux obtained with a null path-length for primaries. 

In the computations presented in figure 5 we have chosen the EBL background of [7], 
nevertheless the flux of secondary protons depends very little on this choice, since the EBL 
effect is restricted to the Lorentz factor range 10^ < L < 2 x 10^ [5]. 

The agreement among the results of SimProp and those of the kinetic approach of [5] 
is not surprising since the former is a direct derivation of the latter in which the photo- 
disintegration process is treated as a fluctuating interaction, through the MC approach dis- 
cussed in section 2. Nevertheless the results presented in this section offer compelling evidence 
of the internal consistency of the computation method presented. 
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Figure 6. Flux at z = in tlie two cases of a pure iron injection with a power law injection 7 = 2.3 
and a maximum acceleration energy E"^""^ = 5 x 10^^ eV (left panel) and of a pure nitrogen injection 
with 7 = 2.0 and E™°-^ = 1.4 x 10^"'^ eV (right panel). The full squares refer to the SimProp all-particle 
spectrum and to the spectrum of secondary nuclei and protons as labeled. The continuous lines refer 
to the corresponding results of [16]. 



Let us conclude this subsection discussing why it is useful to go beyond the kinetic 
approach. The kinetic approach has the important feature of being analytical: fluxes are 
computed mathematically solving a first principles equation [5]. This means that the flux 
of primaries and secondaries is expressed in terms of several integrals that can be computed 
numerically, once the injection spectrum and the sources distribution are specified. In par- 
ticular, the flux of secondary nuclei and nucleons produced in the photo-disintegration chain 
of a primary Aq is obtained by the numerical computation of ^0 nested integrals and this 
computation should be repeated each time the hypothesis on sources (injection and distri- 
bution) are changed. This computation, while it is always feasible numerically, takes some 
time that can be substantially reduced using a MC computation scheme. This follows by the 
fact that, as discussed in section 3, within the SimProp approach it is possible to simulate 
different source distributions and injection spectra without repeating the overall propagation 
of particles. In this sense the MC approach presented here, which is the minimal stochastic 
extension of the kinetic approach, provides a faster computation scheme. Finally, through 
the MC approach of SimProp one takes into account also the intrinsic stochastic nature of 
the photo-disintegration process which is neglected in the kinetic approach. 

4.2 Other MC simulations 

The first code we have chosen for the comparison of SimProp with other MC approaches 
is the one presented in [16]. We used the fluxes reported in ref. [16]. Therefore the same 
injection conditions adopted in [16] have been fixed for SimProp. In particular, two cases of a 
pure injection have been considered: iron nuclei with a power law injection index 7 = 2.3 and 
nitrogen nuclei with 7 = 2.0. The maximum acceleration energy is fixed to E'^°'^ = 5 x 10^^ 
eV in the case of iron and to E"^""^ = 1.4 x 10^^ eV in the case of nitrogen. The results of the 
comparison are shown in figure 6: left panel refers to the case of iron injection and right panel 
to the case of nitrogen. SimProp spectra and the spectra of [16] are normalized to Auger data 
[38] above 10^^'^ eV. From fi gure 6 we can conclude that there is a good agreement among 
the two computations schemes in all spectrum components (primary nuclei, secondary nuclei 
and secondary protons), the largest difference being of the order of 30 % for the flux of 
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intermediate mass nuclei in nitrogen injection at ii^ = 10^^-'^ eV. The results of SimProp and 
by Allard et al. [16] shown in figure 6 are obtained assuming the same model for the EBL 
background given in [7]. 




Figure 7. Flux at z = in the case of a pure iron injection with a power law injection 7 = 2.3 and a 
maximum acceleration energy Eg''"-^ = 5 x 10^^ eV. The full squares refer to the SimProp all-particle 
spectrum and to the spectrum of secondary nuclei and protons as labeled. The continuous lines refer 
to the corresponding results of obtained with CRPropa [28]. 

The comparison between SimProp and CRPropa^ [28] has been done using the same 
choice of injection spectral index and maximum energy as above. The publicly available 
CRPropa framework - which was designed for the simulation of the propagation of nucleons - 
has been recently extended to the propagation of heavy nuclei. The results of the comparison 
are given in figure 7, where the spectra are normalized to Auger data [38] above lO^^'* eV. 
The largest differences (60%) in this comparison are found in the secondary flux of particles 
with 13 < Z < 20 at E = loi8-75 eV. 

The nuclear model adopted in SimProp (see section 3) is simplifled with respect to the 
model used by Allard et al. in [16] and the model used in CRPropa [28]. The good agreement 
of the all-particle spectra demonstrates that a simplifled scheme is effective in producing a 
reliable description of the propagated spectra, especially if we take into account the limited 
mass resolution of the experimental data. 



5 Comparison with the Auger Spectrum 

In this section we compare the spectrum obtained with SimProp with the latest results of 
Auger [38]. This comparison has only illustrative purposes, to show the capabilities of our 
computation scheme. We do not want here to develop a systematic study of the Auger 
observations in terms of spectra, which is outside the scope of this paper. 

To this effect we refer to the spectrum obtained with the simulation code described 
in [16] and calculated for comparison with Auger data presented in the ICRC 2009 [39]: 
therefore we restrict the choice of the injection spectral index to 7 = 2.4 as in this analysis ^. 

In flgure 8 we show the spectra derived for a single component of primaries (Fe) with 
two different values for the iron maximum energy: E^'^^ = 3 x lO^'' eV and E^'^^ = 00. The 

^We thank Simone Riggi, and CRPropa development team, for producing the simulation used in figure 7. 
*For these comparisons, we always generated 400000 events in four rcdshift ranges: 0.0 — 0.02, 0.02 — 
0.2, 0.2 - 1.0, 1.0 - 3.0, using EBL 1. 
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Figure 8. Flux of Fc nuclei injected with a power law with 7 = 2.4 and Eg^"-^' = 3 x 10^° eV for iron 
(left panel) and E^^"-^ — 00 (right panel). The red line refers to the contributions of nuclei (summed 
over iron and all secondaries) , the blue line to the contribution of secondary protons and the grey line 
is the all-particle spectrum. For comparison the Auger combined spectrum [38] is shown with black 
circles. 



sources are assumed to be uniformly distributed in co-moving coordinates, and the all-particle 
spectra are normalized to Auger events above 10^^'^ eV. 

The contributions of nuclei lighter than iron to the all particle spectrum are due to the 
effect of photo-disintegration, that provides at z = a mixture of all secondary protons and 
nuclei with A < 56 in the photo-disintegration chain of iron. As expected, there is remarkable 
difference in the proton fraction at z = depending on the iron maximum acceleration energy. 

From figure 8, we can observe that the different choices of the iron maximum acceleration 
energy have little impact on the all-particle spectra, because of the spectral cut-off due to 
photo-disintegration. However, the proton fraction at high energy is very different in the two 
cases implying some effects of this parameter on the observed elongation rate. 

In figure 9 we plot the fluxes computed with SimProp in the case of 50% injection 
of protons and iron nuclei, and a rigidity dependent energy cut-off (left panel), that is 
E^ax^^^ = ^^«^(Fe)§. 

6 Discussion and Future Development 

In this paper we have presented a new simulation code, SimProp, to simulate the propagation 
of UHE particles in astrophysical backgrounds. The code is based on the analytical scheme of 
[5], modified to take into account possible fluctuations in the photo-disintegration process of 
nuclei. Spectra obtained with SimProp have been successfully checked with spectra obtained 
in the kinetic approach of [5] and with the MC simulation codes of [16] and [28]. 

The approximations presently used in SimProp are mainly related to the nuclear model, 
which is based on a fixed list of nuclei starting from iron down to deuterium. We have 
neglected here the effects due to radioactive decays of nuclei and to the process of photo-pion 
production for nuclei. In future works we will refine SimProp by including also these effects 
and comparing it with other existing simulation codes. Nevertheless, as discussed in section 
4, the good agreement among fluxes obtained with SimProp and with other MC approaches 
shows the tiny dependence of the results on the nuclear model assumed. 

The performances of SimProp depend on the EBL background assumed: choosing a 
simple analytical approximation for this background (parameterization suggested in [7]) the 
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Figure 9. As in figure 8 witli an injection of 50% Fc nuciei and 50% protons, injected witli 7 = 2.4 
and E^''"" = 3 X lO^o eV (left panel), i;™*^^ = 00 (right panel). 



execution time is strongly reduced, as it is shown in table 3 in appendix A. In this case 
SimProp is a very convenient tool to produce enough statistics by obtaining spectra and 
composition observables very quickly. Even using the complete EBL model suggested in 
[7] , the execution time of SimProp is reasonable enabling a fast analysis of the results (see 
appendix A). 

Future developments of SimProp are planned. As a first enhancement, the introduction 
of the photo-pion production process for nuclei will be taken into account. Three-dimensional 
effects caused by the granularity of the actual source distribution and the effects of magnetic 
fields in the propagation of nuclei are not presently included. The results presented here 
are all obtained assuming a uniform distribution of sources in co-moving coordinates, a case 
which gives a flux independent of the magnetic field [40]. We will also improve SimProp to 
carry on a systematic study of the effects of a sparse distribution of sources and of galactic 
as well as (possible) extra-galactic magnetic fields. 

Finally, in the present version of SimProp the treatment of secondary photons and 
neutrinos produced by the propagation of protons and nuclei is not implemented. It has to 
be noticed, however, that the main ingredients to determine the fluxes of such secondaries 
are already present in the outputs of our code, being just the production energy and red-shift 
of protons. Therefore we will soon consider also this extension of our simulation code. 

A Performances 

In this section the averaged execution times per event (seconds) are reported for different 
cases. The tests have been performed using a virtual machine QEMU Virtual CPU at 
2.27GHz with 2 GB RAM. The reults are shown in the following tables. 

The main effect on the performances of SimProp is given by the energy of the simulated 
events while the influence of the redshift range is scarcely affecting the computing times (see 
table 1). 

In table 2 it is shown the mean execution times for different primary masses at injection. 
The redshift range does not influence too much the performances of SimProp, while there is 
an increase of the computation time with the primary nuclear type. 

The last test was performed to compare the performance of SimProp depending on 
the parametrization of the background. Two different parameterization have been used and 
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described in section 2. The comparison in table 3 shows an increase of the computation time 
of about 40 times going from EBL 2 to EBL 1 . Also for the simulation of background EBL 
1 the simulation time weakly depends on the redshift range. 





< z < 0.2 


0.2 < z < 1 


1 < z < 3 


16 < logio(^/eV) < 18 


0.007 


0.009 


0.007 


18 < logio(^/eV) < 20 


0.38 


0.95 


0.96 


20 < logio(^/eV) < 21.5 


2.03 


1.83 


1.48 



Table 1. Execution times (seconds) per event (Fe injection) for different energy ranges of the primary 
and redshift ranges of the distribution of sources. 





< z < 0.2 


0.2 < z < 1 


1< z <3 


A = 4 


0.18 


0.19 


0.19 


A = U 


0.38 


0.39 


0.37 


^ = 56 


1.1 


1.3 


1.1 



Table 2. Execution times (seconds) per event for different species at injection corresponding to 
an energy range of the primary 17 < log]^o(i?/eV) < 22.5 and to different redshift ranges of the 
distribution of sources. 





< z < 0.2 


0.2 < z < 1 


1< z < 3 


EBL 1 


41.5 


44.6 


44.5 


EBL 2 


1.1 


1.3 


1.1 



Table 3. Execution times (seconds) per event for different paramctrizations of the EBL back- 
ground distribution of photons (see text) corresponding to an energy range of the primary Fe 
17 < log2o(-E/eV) < 22.5 and to different redshift ranges of the distribution of sources. 
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